#pragma once

#include "Matrix.h"
#include "RBT.h"
#include <deque>
#include <string>

/**
 * @file Mesh.h
 * @author xingyifan
 * @date 2023-01-14 14:45
 *
 * @description: 异常类型
 */

#define DEBUG_NULL_VERTEX 0		// 空顶点
#define DEBUG_NULL_HALFEDGE 1	// 空半边
#define DEBUG_NULL_FACE 2		// 空面
#define DEBUG_EXIST_TWIN 3		// 已经存在配对
#define DEBUG_ERROR_EDGE_MESH 4 // 网格半边错误

/**
 * @file Mesh.h
 * @author xingyifan
 * @date 2023-01-14 14:59
 *
 * @description: 异常描述
 */

#define DEBUG_NULL_VERTEX_DETAIL "Some vertex is nullptr.\n"
#define DEBUG_NULL_HALFEDGE_DETAIL "Some half-edge is nullptr.\n"
#define DEBUG_NULL_FACE_DETAIL "Some face is nullptr.\n"
#define DEBUG_EXIST_TWIN_DETAIL "Some twin exists.\n"
#define DEBUG_ERROR_EDGE_MESH_DETAIL "Some edge has more than two half-edge.\n"

namespace CG
{
	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2023-01-04 11:58
	 *
	 * @description: 基础数据类型
	 */

	enum class CGElement
	{
		vertex,
		half_edge,
		face
	};

	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2023-01-04 11:59
	 *
	 * @description: 包围盒
	 */

	struct bounding_box
	{
		double left, right;
		double top, bottom;
		double zNear, zFar;
	};

	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2022-12-15 13:22
	 *
	 * @description: 三维点/向量
	 */

	class vector3f
	{
	public:
		double x, y, z;

	public:
		vector3f(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}

		// 基本运算
		vector3f operator*(double t) const { return {x * t, y * t, z * t}; }
		vector3f operator/(double t) const
		{
			DEBUG_ASSERT(fabs(t) > 1e-16, DEBUG_ZERO_DIVISION, vector3f())
			return {x / t, y / t, z / t};
		}
		vector3f operator-(const vector3f &v) const { return {x - v.x, y - v.y, z - v.z}; }
		vector3f operator+(const vector3f &v) const { return {x + v.x, y + v.y, z + v.z}; }
		friend vector3f operator-(const vector3f &v) { return {-v.x, -v.y, -v.z}; }
	};

	class point3f
	{
	public:
		double x, y, z;

	public:
		point3f(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}

		vector3f operator-(const point3f &p) const { return {x - p.x, y - p.y, z - p.z}; }							// 顶点相减得到向量
		friend point3f operator+(const point3f &p, const vector3f &v) { return {p.x + v.x, p.y + v.y, p.z + v.z}; } // 顶点 + 向量 = 顶点
		friend point3f operator+(const vector3f &v, const point3f &p) { return {p.x + v.x, p.y + v.y, p.z + v.z}; } // 向量 + 顶点 = 顶点
	};

	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2023-01-04 12:03
	 *
	 * @description: 4 维齐次坐标
	 */

	struct coord4f
	{
		double x, y, z, w;
	};

	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2023-01-04 12:07
	 *
	 * @description: RGB 空间
	 */

	class rgb_space
	{
	public:
		double r, g, b;

	public:
		rgb_space(double r = 0, double g = 0, double b = 0) : r(r), g(g), b(b) {}

		// 颜色直接叠加
		rgb_space operator+(const rgb_space &c1) const
		{
			rgb_space c = {c1.r + r, c1.g + g, c1.b + b};
			c.r = (c.r > 1) ? 1 : c.r;
			c.g = (c.g > 1) ? 1 : c.g;
			c.b = (c.b > 1) ? 1 : c.b;
			return c;
		}
	};

	// 材质结构
	struct material
	{
		double alpha;	// 镜面反射指数
		double R[3][3]; // 反射率矩阵 a d s ：环境光、漫反射、镜面反射
	};

	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2023-01-04 12:12
	 *
	 * @description: 基本数据操作
	 */

	// 颜色线性插值
	inline rgb_space linear(rgb_space c1, rgb_space c2, double t)
	{
		return {c1.r * (1 - t) + c2.r * t,
				c1.g * (1 - t) + c2.g * t,
				c1.b * (1 - t) + c2.b * t};
	}

	// 点线性插值
	inline point3f linear(point3f p1, point3f p2, double t, double a = 0, double b = 1)
	{
		return {(p1.x * (b - t) + p2.x * (t - a)) / (b - a),
				(p1.y * (b - t) + p2.y * (t - a)) / (b - a),
				(p1.z * (b - t) + p2.z * (t - a)) / (b - a)};
	}

	// 计算点的距离
	inline double distance(point3f p1, point3f p2)
	{
		return sqrt(pow(p1.x - p2.x, 2) + pow(p1.y - p2.y, 2) + pow(p1.z - p2.z, 2));
	}

	// 向量点积
	inline double dot(vector3f v1, vector3f v2)
	{
		return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
	}

	// 向量叉积
	inline vector3f cross(vector3f v1, vector3f v2)
	{
		return {v1.y * v2.z - v1.z * v2.y,
				-v1.x * v2.z + v1.z * v2.x,
				v1.x * v2.y - v1.y * v2.x};
	}

	// 行列式 (v1 x v2) * v3
	inline double det(vector3f v1, vector3f v2, vector3f v3)
	{
		return dot(cross(v1, v2), v3);
	}

	// 向量长度
	inline double length(vector3f p)
	{
		return sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
	}

	// 向量规范化
	inline vector3f normalize(vector3f p)
	{
		double len = length(p);
		DEBUG_ASSERT(len > 1e-16, DEBUG_ZERO_DIVISION, vector3f())
		return p / len;
	}

	// 判断三点是否共线
	inline bool onLine(point3f p1, point3f p2, point3f p3)
	{
		vector3f vec1 = p1 - p2;
		vector3f vec2 = p1 - p3;

		// 叉积足够大，则认为不共线
		if (length(cross(vec1, vec2)) > 1e-16)
			return false;
		return true;
	}

	// 判断四点是否共面
	inline bool onFace(point3f p1, point3f p2, point3f p3, point3f p4)
	{
		vector3f vec1 = p1 - p2;
		vector3f vec2 = p1 - p3;
		vector3f vec3 = p1 - p4;

		// 计算三个向量的行列式，如果足够大，则认为不共面
		if (fabs(det(vec1, vec2, vec3)) > 1e-16)
			return false;
		return true;
	}

}; // namespace CG

namespace CG
{
	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2023-01-04 15:54
	 *
	 * @description: 四元数类型
	 */

	class quaternion
	{
	public:
		double s;
		vector3f v;

	public:
		// 常规初始化，分量初始化、默认初始化，纯四元数初始化
		quaternion(double a, double b, double c, double d) : s(a), v({b, c, d}) {}
		quaternion(double s = 0, vector3f v = {0, 0, 0}) : s(s), v(v) {}
		quaternion(vector3f v) : s(0), v(v) {}

		// 获得绕轴 v 旋转 theta 角的四元数
		quaternion(vector3f v, double theta) : s(cos(theta / 2)), v(normalize(v) * sin(theta / 2)) {}

		// 四元数基本运算
		quaternion operator*(double t) const { return quaternion(s * t, v * t); }
		quaternion operator/(double t) const
		{
			DEBUG_ASSERT(fabs(t) > 1e-16, DEBUG_ZERO_DIVISION, quaternion())
			return quaternion(s / t, v / t);
		}
		quaternion operator+(const quaternion &q) const { return quaternion(s + q.s, v + q.v); }
		quaternion operator-(const quaternion &q) const { return quaternion(s - q.s, v - q.v); }
		quaternion operator*(const quaternion &q) const { return quaternion(s * q.s - dot(v, q.v), q.v * s + v * q.s + cross(v, q.v)); }
		quaternion operator/(const quaternion &q) const { return (*this) * quaternion(q.s, -q.v) / pow(sqrt(q.s * q.s + pow(length(q.v), 2)), 2); }
	};

	// 四元数长度
	inline double length(quaternion q)
	{
		return sqrt(q.s * q.s + pow(length(q.v), 2));
	}

	// 四元数规范化
	inline quaternion normalize(quaternion q)
	{
		int len = length(q);
		DEBUG_ASSERT(len, DEBUG_ZERO_DIVISION, quaternion())
		return q / len;
	}

	// 获得从轴 a 转动到轴 b 的四元数
	inline quaternion rotate(vector3f a, vector3f b)
	{
		// 确保两个轴长度不为零
		DEBUG_ASSERT(length(a) > 1e-16 && length(b) > 1e-16, DEBUG_ZERO_DIVISION, quaternion())

		// 先把 a 转动到 z 轴：先绕 z 轴转，然后绕 y 轴转
		quaternion q1(1), q2(1);
		double proj = sqrt(pow(a.x, 2) + pow(a.y, 2));

		// 如果已经在 z 轴，就不动；否则计算两个四元数
		if (proj > 1e-16)
		{
			q1 = quaternion({0, 0, 1}, -acos(a.x / proj));
			q2 = quaternion({0, 1, 0}, -acos(a.z / length(a)));
		}

		// 然后把 b 同样转动到 z 轴
		quaternion q3(1), q4(1);
		proj = sqrt(pow(b.x, 2) + pow(b.y, 2));

		// 如果已经在 z 轴，就不动；否则计算两个四元数
		if (proj > 1e-16)
		{
			q3 = quaternion({0, 0, 1}, -acos(b.x / proj));
			q4 = quaternion({0, 1, 0}, -acos(b.z / length(b)));
		}
		return quaternion(1) / q3 / q4 * q2 * q1;
	}
}; // namespace CG

namespace CG
{
	// 网格基本结构
	struct vertex3f;
	struct half_edge;
	struct face;

	// 网格类型
	enum class MeshType
	{
		empty,
		triangular,
		rectangular,
		mixed
	};

	// 3 维顶点
	struct vertex3f
	{
		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-11-25 20:57
		 *
		 * @description: 基础信息
		 */

		point3f pos;  // 坐标
		half_edge *e; // 以它为起点的一条半边

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-11-25 20:58
		 *
		 * @description: 附加信息
		 */

		int identity; // 顶点 id 标记
	};

	// 面片
	struct face
	{
		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-11-25 20:58
		 *
		 * @description: 基础信息
		 */

		half_edge *e; // 面边界上的一条半边

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-11-25 20:58
		 *
		 * @description: 附加信息
		 */

		int identity; // 面标记
	};

	// 半边结构
	struct half_edge
	{
		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-12-17 22:25
		 *
		 * @description: 基础信息
		 */

		vertex3f *v;	 // 起点
		half_edge *twin; // 反向边，不一定存在
		face *incFace;	 // 沿着半边方向左边的面，这个面一定存在
		half_edge *prev; // 上一条半边，一定存在，在面上
		half_edge *next; // 下一条半边，一定存在，在面上

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-12-17 22:25
		 *
		 * @description: 附加信息
		 */

		int identity; // 半边标记
	};

	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2023-01-10 16:58
	 *
	 * @description: 补充算法类
	 */

	class valid_pair
	{
	public:
		half_edge *e; // 保存半边
		LA::Vector v; // 记录最优解，前三位是坐标，最后一位存放误差

	public:
		// RBT 需要默认构造函数创建 nil 节点
		valid_pair() : e(nullptr), v(LA::Vector(4)) {}
		valid_pair(half_edge *e, LA::Vector v) : e(e), v(v) {}

		// 根据误差值比较，用于红黑树插入
		bool operator<(const valid_pair &vp) const { return v[3] < vp.v[3]; }
	};

	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2022-12-04 11:38
	 *
	 * @description: 基于半边结构的网格
	 */

	class Mesh
	{
	public:
		Mesh() : m_type(MeshType::triangular)
		{
			// 需要默认构造材质，默认为纯白色
			m_material.alpha = 100;
			for (int i = 0; i < 3; i++)
				for (int j = 0; j < 3; j++)
					m_material.R[i][j] = 1;

			// 默认包围盒为空
			GetBox();
		}
		// 最后释放内存
		~Mesh()
		{
			for (auto it = m_vList.begin(); it != m_vList.end(); ++it)
				if ((*it) != nullptr)
					delete *it;
			for (auto it = m_eList.begin(); it != m_eList.end(); ++it)
				if ((*it) != nullptr)
					delete *it;
			for (auto it = m_fList.begin(); it != m_fList.end(); ++it)
				if ((*it) != nullptr)
					delete *it;
		}

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-11-25 13:16
		 *
		 * @description: 初始化构造操作
		 */

		vertex3f *CreateVertex(point3f p, bool testSame = false);										// 创建顶点，并返回指针（可选检查顶点是否重合，如果重合，将返回重合顶点）
		face *CreateFace(vertex3f *v0, vertex3f *v1, vertex3f *v2, bool autoTwin = true);				// 创建三角面并返回，默认自动判定衔接相邻面
		face *CreateFace(vertex3f *v0, vertex3f *v1, vertex3f *v2, vertex3f *v3, bool autoTwin = true); // 创建四边形面，默认自动判定衔接相邻面
		face *CreateFace(vertex3f **vList, int n, bool autoTwin = true);								// 创建多边形面
		void CreateTwin(face *f);																		// 为一个面创建连接，会选择面的顶点相邻的面尝试建立连接

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-11-25 20:34
		 *
		 * @description: 获取基本网格信息
		 */

		size_t VertexNum() const { return m_vList.size(); }
		size_t FaceNum() const { return m_fList.size(); }
		size_t HalfEdgeNum() const { return m_eList.size(); }

		bounding_box GetBox(); // 获得网格的包围盒，同时会储存获得的包围盒

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-11-25 13:16
		 *
		 * @description: 结构删除操作
		 */

		bool DeleteVertex(vertex3f *v);		 // 删除顶点和相邻的面及其半边
		bool DeleteHalfEdge(half_edge *e);	 // 删除半边和相邻的面，将目标顶点合并到起点，目标顶点将会孤立
		bool DeleteFace(face *f);			 // 删除面和对应的半边
		bool DeleteIsolated(CGElement type); // 清除孤立元素

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-12-15 14:07
		 *
		 * @description: 附加判断
		 */

		bool HasVertexNormal();	  // 是否有顶点法向
		bool IsConvexHull();	  // 判断网格是否是凸包
		bool IsInMesh(point3f p); // 判断点是否在网格包围中

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-12-15 13:27
		 *
		 * @description: 网格结构方法
		 */

		bool StitchFace(face *f1, face *f2);					  // “缝合”两个面，返回是否出错
		bool OnFace(point3f p, face *f);						  // 判断一个点是否在面上
		bool AboveFace(point3f p, face *f);						  // 判断一个点是否在面的上方
		vector3f Gouraud(vertex3f *v);							  // Gouraud 插值：计算顶点的平均法向
		vector3f FaceNormal(face *f);							  // 计算单位面法向
		vertex3f *Destination(half_edge *e);					  // 获得边的目标点
		std::deque<vertex3f *> VertexOnFace(face *f);			  // 获得一个面上的所有顶点
		std::deque<half_edge *> HalfEdgeOnFace(face *f);		  // 获得一个面上的所有半边
		std::deque<half_edge *> CommonOnFace(face *f0, face *f1); // 获取两个面的公共半边
		std::deque<vertex3f *> AdjacentVertex(vertex3f *v);		  // 获得 v 相邻的所有顶点
		std::deque<half_edge *> AdjacentHalfEdge(vertex3f *v);	  // 获得以 v 为起点的半边
		std::deque<face *> AdjacentFace(vertex3f *v);			  // 获得 v 相邻的所有面
		std::deque<face *> AdjacentFace(face *f);				  // 获得 f 相邻的所有面
		std::deque<vertex3f *> MeshProject(vector3f pv);		  // 取网格点沿着 pv 方向的投影多边形顶点
		std::deque<face *> TriangulateFace(face *f);			  // 将网格面三角化，返回三角化后的面列表
		LA::Matrix QEMatrix(vertex3f *v);						  // 获得顶点的二次误差矩阵
		LA::Vector MinVector(half_edge *e);						  // 获得边上顶点收缩的最优解
		bool Simplify(int times);								  // 曲面简化算法 QEM，不考虑非流形网格：有效对总是边。能否简化成功

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2023-01-10 20:56
		 *
		 * @description: 网格属性设置
		 */

		void SetMaterial(double M[][3], double alpha); // 设置表面系数：反射系数和镜面反射系数

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-12-25 19:33
		 *
		 * @description: 保存和读取操作
		 */

		bool SaveMesh(const std::string &filename);						   // 保存网格顶点、面和邻面关系
		bool LoadMesh(const std::string &filename, bool testSame = false); // 加载网格，可选是否排除重合顶点

	public:
		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-12-15 13:55
		 *
		 * @description: 基本数据结构
		 */

		std::deque<vertex3f *> m_vList;	 // 顶点列表
		std::deque<half_edge *> m_eList; // 半边列表
		std::deque<face *> m_fList;		 // 面列表
		std::deque<vector3f> m_nList;	 // 法向列表

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2023-01-10 20:51
		 *
		 * @description: 进一步的网格属性
		 */

		std::string m_name;	 // 网格名
		MeshType m_type;	 // 网格类型
		material m_material; // 网格材质
		bounding_box m_box;	 // 网格包围盒

		/**
		 * @file Mesh.h
		 * @author xingyifan
		 * @date 2022-12-15 13:55
		 *
		 * @description: 附加信息，这部分内容只是存放在 Mesh 中，并不作为基本属性使用
		 */

		std::deque<rgb_space> m_cList;	// 颜色列表
		std::deque<LA::Matrix> m_QList; // 误差矩阵列表
	};

	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2022-12-17 22:19
	 *
	 * @description: 增量凸包算法
	 */

	bool tetrahedron(Mesh *mesh);								  // 尝试用网格点创建四面体，返回是否成功
	bool convexHull(Mesh *mesh);								  // 取网格点的凸包，返回是否成功
	bool isInHull(point3f p, std::deque<point3f> &hull);		  // 检查点是否在凸包内部
	std::deque<point3f> convexHull(std::deque<point3f> &polygon); // 计算并返回平面多边形的凸包

	/**
	 * @file Mesh.h
	 * @author xingyifan
	 * @date 2022-12-26 11:41
	 *
	 * @description: 常见曲面生成
	 */

	bool sphere(Mesh *mesh, point3f center, vector3f axis, double a, double b, double c, int div = 20);			// 球面：中心位置、中轴、三个半长轴，分割份数
	bool cylinder(Mesh *mesh, point3f center, vector3f axis, double r, double h, int divr = 40, int divh = 40); // 柱面：中心位置，中轴，半径，高度，分割份数
	bool cone(Mesh *mesh, point3f origin, vector3f axis, double theta, double h, int divr = 40, int divh = 40); // 锥面：中心位置，中轴，半角，高度，分割份数
	bool ring(Mesh *mesh, point3f center, vector3f axis, double r, double R, int divR = 40, int divr = 40);		// 圆环面：中心位置，垂直轴，内径，外径，圆的份数，每个圆有多少个点

}; // namespace CG
